source('/home/jc148322/scripts/libraries/cool_functions.r')
library(SDMTools)
climate.dir='/home/jc165798/Climate/CIAS/AWT/250m/baseline.76to05/bioclim/'
wd='/home/jc148322/frog_declines/SDM/';setwd(wd)
out.dir='/home/jc148322/frog_declines/SDM/boxplots/'; dir.create(out.dir)

timeline=c('Pre_decline','Post_decline')

clim.vars = paste('bioclim_',sprintf('%02i',c(1,4,5,6,12,15,16,17)),sep='')

species=list.files(paste(wd,timeline[2],'/bioclim_1_4_5_6_12_15_16_17/models/',sep=''))

base.asc=read.asc.gz('/home/jc148322/frog_declines/SDM/Pre_decline/bioclim_12_15/models/acuti/output/ascii/current.76to05.asc.gz')		

pos=make.pos(base.asc)
clim.vars = paste('bioclim_',sprintf('%02i',c(1,4,5,6,12,15,16,17)),sep='')
for (clim.var in clim.vars) { cat(clim.var, '\n')
	pos[,clim.var]=extract.data(cbind(pos$lon,pos$lat),read.asc.gz(paste(climate.dir,clim.var,'.asc.gz',sep='')))
}

for (spp in species) {cat(spp,'\n')
	for (tt in timeline) { cat(tt,'\n')
		md=paste(wd, tt, '/bioclim_1_4_5_6_12_15_16_17/models/',spp,'/output/',sep='');setwd(md)
		
		threshold=read.csv('maxentResults.csv',as.is=TRUE)
		threshold=threshold$Balance.training.omission..predicted.area.and.threshold.value.logistic.threshold
		
		tasc=read.asc.gz(paste('ascii/current.76to05.asc.gz', sep=''))
		tasc[which(tasc<=threshold)]=0
		
		spdata=pos
		spdata$suitability=extract.data(cbind(spdata$lon,spdata$lat),tasc)
		spdata$suitability[which(spdata$suitability>0)]=1
		spdata[which(spdata$suitability==0),grep('bioclim',colnames(spdata))]=NA
		spdata=na.omit(spdata[,grep('bioclim',colnames(spdata))])
		# colnames(spdata)=paste(tt,'_',sprintf(colnames(spdata)),sep='')

		assign(tt,spdata)
	}

	setwd(out.dir)
	for (clim.var in clim.vars) { cat(clim.var,'\n')
		tdata=list(Pre_decline[,clim.var],Post_decline[,clim.var])
		png(paste(out.dir,spp,'_',clim.var,'.png',sep=''))
			boxplot(tdata)
		dev.off()

	}	
}	

out.dir='/home/jc148322/frog_declines/SDM/boxplots_occur/'; dir.create(out.dir)
clim.vars=paste('bioclim_',sprintf('%02i',c(1:19)),sep='')
for (clim.var in clim.vars){ cat(clim.var,'\n')
	tasc=read.asc.gz(paste(climate.dir,clim.var,'.asc.gz',sep=''))
	assign(clim.var,tasc)
}

for (spp in species) {cat(spp,'\n')
	for (tt in timeline) { cat(tt,'\n')
		md=paste(wd, tt, '/bioclim_1_4_5_6_12_15_16_17/models/',spp,'/',sep='');setwd(md)
		
		occur=read.csv('occur.csv',as.is=TRUE)
		occur=occur[1:3]
		for (clim.var in clim.vars) { cat(clim.var,'\n') #cycle through each of hte climate variables
			occur[,clim.var] = extract.data(cbind(occur$long,occur$lat), get(clim.var)) #append the climate data
		}
		assign(tt,occur)
	}

	setwd(out.dir)
	for (clim.var in clim.vars) { cat(clim.var,'\n')
		tdata=list(Pre_decline[,clim.var],Post_decline[,clim.var])
		png(paste(out.dir,spp,'_',clim.var,'.png',sep=''))
			boxplot(tdata)
		dev.off()

	}	
}	

	